<!--
Plotting demo for cubic spline interpolation.

To modify the demo, adjust the commented-out lines in the code below, or
customize in any other way.

Eli Bendersky [https://eli.thegreenplace.net]
This code is in the public domain.
-->
<!DOCTYPE html>
<svg id="chart" width="500" height="300"></svg>

<script type="module">
    import * as d3 from "https://cdn.jsdelivr.net/npm/d3@7/+esm";
    import { buildSplineEquations } from "./spline.js";
    import { solve } from "./eqsolve.js";

    // origXs, origYs is the original set of points we're going to be
    // interpolating between. They can be any length, as long as
    // origXs.length == origYs.length

    // const origXs = [0, 1, 2];
    // const origYs = [1, 3, 2];

    // const origXs = [0, 1, 2, 3, 4];
    // const origYs = [21, 24, 24, 18, 16];

    // sinc function sampled on a few uniformly-spaced points.
    const origXs = linspace(-10, 10, 7);
    const origYs = origXs.map(calcSinc);

    // Interpolate 200 points between the original points, using cubic spline
    // interpolation.
    let [pxs, pys] = doInterpolate(origXs, origYs, 200);

    // Plotting code with D3.
    const xMin = d3.min([...origXs]), yMin = d3.min([...origYs]);
    const xMax = d3.max([...origXs]), yMax = d3.max([...origYs]);

    const margin = { top: 20, right: 20, bottom: 40, left: 40 };
    const width = document.querySelector('#chart').getAttribute('width') - margin.left - margin.right;
    const height = document.querySelector('#chart').getAttribute('height') - margin.top - margin.bottom;

    const svg = d3.select("#chart")
        .append("g")
        .attr("transform", `translate(${margin.left}, ${margin.top})`);

    const xScale = d3.scaleLinear().domain([xMin - 1, xMax + 1]).range([0, width]);
    const yScale = d3.scaleLinear().domain([yMin, yMax]).range([height, 0]);

    const xAxis = d3.axisBottom(xScale);
    const yAxis = d3.axisLeft(yScale);

    svg.append("g")
        .attr("transform", `translate(0, ${height})`)
        .call(xAxis);

    svg.append("g")
        .call(yAxis);

    // Linear connections
    const linear = d3.line()
        .x((d, i) => xScale(origXs[i]))
        .y(d => yScale(d));
    svg.append("path")
        .attr("d", linear(origYs))
        .attr("stroke", "gray")
        .attr("stroke-width", "1")
        .attr("fill", "none");

    // Interpolated line
    const interp = d3.line()
        .x((d, i) => xScale(pxs[i]))
        .y(d => yScale(d));
    svg.append("path")
        .attr("d", interp(pys))
        .attr("stroke", "blue")
        .attr("stroke-width", "2")
        .attr("fill", "none");

    // Dots on the original points
    svg.selectAll(".dot")
        .data(origYs)
        .enter().append("circle")
        .attr("stroke", "red")
        .attr("fill", "red")
        .attr("class", "dot")
        .attr("cx", (d, i) => xScale(origXs[i]))
        .attr("cy", d => yScale(d))
        .attr("r", 3);

    // doInterpolate uses cubic spline interpolation to create N new points
    // between xs and calculates their ys, returning [pxs, pys] - the (x,y)
    // coords of the interpolated points.
    function doInterpolate(xs, ys, N) {
        // Perform interpolation on xs, ys to get the coefficients of the splines.
        let [A, b] = buildSplineEquations(xs, ys);
        let coeffs = solve(A, b);
        console.log(coeffs);

        // Create N points linearly spaced between the min and max of xs, and
        // calculate the corresponding py for each px using the appropriate curve.
        let pxs = linspace(Math.min(...xs), Math.max(...xs), N);

        let pys = Array(N).fill(0);
        for (let i = 0; i < N; i++) {
            let px = pxs[i];
            // Find the number of the curve for px, based on which points from
            // xs it's between. Can be done more efficiently with binary
            // search, but this is good enough for a demo.
            let curveIndex = -1;
            for (let j = 0; j < xs.length - 1; j++) {
                // is px between xs[j] and xs[j+1]? If yes, we found the curve!
                if (px >= xs[j] && px <= xs[j + 1]) {
                    curveIndex = j;
                    break;
                }
            }
            if (curveIndex < 0) {
                alert(`curve index not found for xs[${i}]=${xs[i]}`);
            }

            // With the curve index in hand, we can calculate py based on the
            // relevant curve coefficients from coeffs.
            let [a, b, c, d] = coeffs.slice(curveIndex * 4, curveIndex * 4 + 4);
            pys[i] = a * px ** 3 + b * px ** 2 + c * px + d;
        }

        return [pxs, pys];
    }

    // linspace returns an array of numPoints values distributed linearly in
    // the (inclusive) range [start,end], just like Numpy's linspace.
    function linspace(start, end, numPoints) {
        if (numPoints === undefined || numPoints < 2) {
            return [start, end];
        }

        const step = (end - start) / (numPoints - 1);
        return new Array(numPoints).fill(null).map((_, i) => start + i * step);
    }

    // calcSinc calculates the sinc(x) function and returns a y value.
    // https://en.wikipedia.org/wiki/Sinc_function
    function calcSinc(x) {
        if (x == 0) {
            return 1
        } else {
            return Math.sin(Math.PI * x) / (Math.PI * x);
        }
    }
</script>